PUSHING  MEASUREMENT  TO  THE  ULTIMATE  STOCHASTIC 
LIMIT:  THE  STOCHASTIC  DYNAMICS  OE  ELUID-COUPLED 

NANOCANTILEVERS 
AEOSR  EA9550-07- 1-0222 
EINAL  REPORT 
Mark  Paul 

Department  of  Mechanical  Engineering 
Virginia  Tech 


Abstract 

We  have  developed  a  fundamental  understanding  of  nanoseale  fluid  dynamies  for  fluid- 
based  teehnologies  with  unpreeedented  eapabilities.  Using  analyties  and  numeries  we 
have  investigated  the  Brownian  driven,  and  externally  driven,  dynamies  of  miero  and 
nanoseale  elastie  objeets  (sueh  as  eantilevers  and  beams)  in  a  viseous  fluid  over  a  wide 
range  of  system  parameters  and  for  a  number  of  experimentally  important  eonfigurations. 
We  developed  an  approaeh  to  eompute  the  Brownian  or  externally  driven  dynamics  using 
a  single  deterministic  computation  that  can  be  performed  on  a  personal  workstation.  Ther¬ 
mal  motion  is  computed  using  the  fluctuation-dissipation  theorem  and  externally  driven 
dynamics  using  transfer  function  theory.  We  quantify  the  effects  of  the  cantilever  and 
beam  geometry  upon  their  dynamics,  the  role  of  nearby  bounding  surfaces,  the  increased 
frequency  and  quality  factors  when  using  the  higher  flexural  modes,  and  build  a  physical 
understanding  of  the  fluid  correlated  motion  of  an  array  of  elastic  objects. 


1  Computing  the  Stochastic  and  Driven  Dynamics  of  Elas¬ 
tic  Objects  in  Fluid 

The  stochastic  dynamics  of  micron  and  nanoseale  cantilevers  driven  by  thermal  or  Brown¬ 
ian  motion  can  be  quantified  using  strictly  deterministic  calculations.  This  is  accomplished 
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using  the  fluctuation-dissipation  theorem  since  the  cantilever  remains  near  thermodynamic 
equilibrium  [1,2].  As  part  of  this  research  project  we  have  extended  this  approach  to  the 
experimentally  important  case  of  determining  the  stochastic  dynamics  of  the  angle  of  the 
cantilever  tip.  The  following  discussion  has  been  reported  in  Refs.  [3,4]. 


The  autocorrelation  of  equilibrium  fluctuations  in  cantilever  displacement  can  be  deter¬ 
mined  from  the  deterministic  response  of  the  cantilever  to  the  removal  of  a  step  force 
from  the  tip  of  the  cantilever  (i.e.  a  transverse  point  force  removed  from  the  distal  end  of 
the  cantilever).  If  this  force  f{t)  is  given  by 
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Fq  for  t  <  0 
0  for  f  >  0, 

where  t  is  time  and  Fq  is  the  magnitude  of  the  force,  then  the  autocorrelation  of  the  equi¬ 
librium  fluctuations  in  the  displacement  of  the  cantilever  tip  is  given  directly  by 
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where  is  Boltzmann’s  constant,  T  is  the  temperature,  and  ()  is  an  equilibrium  ensemble 
average.  In  our  notation  lower  case  letters  represent  stochastic  variables  (ui{t)  is  the 
stochastic  displacement  of  the  cantilever  tip)  and  upper  case  letters  represent  deterministic 
variables  represents  the  deterministic  ring  down  of  the  cantilever  tip  due  to  the  step 
force  removal).  The  spectral  properties  of  the  stochastic  dynamics  are  given  by  the  Fourier 
transform  of  the  autocorrelation. 


The  thermodynamic  approach  is  valid  for  any  conjugate  pair  of  variables  [5] .  For  example, 
it  is  common  in  experiment  to  use  optical  techniques  to  measure  the  angle  of  the  cantilever 
tip  as  a  function  of  time  [6].  It  has  also  been  proposed  to  use  piezoresistive  techniques  to 
measure  voltage  as  a  function  of  time  [7].  The  thermodynamic  approach  remains  valid  for 
these  situations  by  choosing  the  correct  conjugate  pair  of  variables. 


We  have  also  explored  the  stochastic  dynamics  of  the  angle  of  the  cantilever  tip.  In  this 
case,  the  angle  of  the  cantilever  tip  is  conjugate  to  a  step  point- torque  applied  to  the  can¬ 
tilever  tip.  If  this  torque  is  given  by 


r(t)  = 


(3) 


To  for  t  <  0 
0  for  t  >  0, 

where  Tq  is  the  magnitude  of  the  step  torque,  then  the  autocorrelation  of  equilibrium  flue 
tuations  in  cantilever  tip-angle  6{t)  is  given  by 
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L(/im)  w(/im)  /i(/im)  k  (N/m)  kt  (N-m/rad)  /o  (kHz) 


(1)  197  29  2  1.3  1.6  x  10"®  71 

(2)  140  15.6  0.6  0.1  8.9  x  10-i°  38 

Table  1:  Summary  of  the  cantilever  geometries  and  material  properties.  (1)  The  rect¬ 
angular  cantilever.  (2)  The  V-shaped  cantilever  used  is  the  commercially  available  Veeco 
MLCT  Type  E  microlever  that  is  used  in  AFM  [8].  The  geometry  is  given  by  the  cantilever 
length  L,  width  w,  and  height  h.  For  the  V-shaped  cantilever  the  total  length  between  the 
two  arms  at  the  base  is  b  =  161.64/rm.  The  cantilever  spring  constant  k,  torsional  spring 
constant  kt,  and  resonant  frequency  in  vacuum  /o  are  determined  using  finite  element  nu¬ 
merical  simulations.  The  cantilevers  are  immersed  in  water  with  density  pi  =  997  kg/m® 
and  dynamic  viscosity  p  =  8.59  x  10“^  kg/m-s. 

Here  0i(t)  represents  the  deterministic  ring  down,  as  measured  by  the  tip-angle,  resulting 
from  the  removal  of  a  step  point-torque.  Again,  the  Fourier  transform  of  the  autocorrela¬ 
tion  yields  the  noise  spectrum. 

A  powerful  aspect  of  this  approach  is  that  it  is  possible  to  use  deterministic  numerical  sim¬ 
ulations  to  determine  Ui{t)  and  0i(t)  for  the  precise  cantilever  geometries  and  conditions 
of  experiment.  This  includes  the  full  three-dimensionality  of  the  dynamics  which  are  not 
accounted  for  in  available  theoretical  descriptions.  The  numerical  results  can  be  used  to 
guide  the  development  of  more  accurate  theoretical  models. 

1.1  The  stochastic  dynamics  of  cantilever  tip-deflection  and  tip-angle 

The  stochastic  dynamics  of  the  cantilever  tip-displacement  Ui{t)  and  that  of  the  tip-angle 
9i{t)  yield  interesting  differences.  Using  the  thermodynamic  approach,  insight  into  these 
differences  can  be  gained  by  performing  a  mode  expansion  of  the  cantilever  using  the 
initial  deflection  required  by  the  deterministic  calculation.  The  two  cases  of  a  tip-force 
and  a  tip-torque  result  in  a  significant  difference  in  the  mode  expansion  coefficients  which 
can  be  directly  related  to  the  resulting  stochastic  dynamics. 

For  small  deflections  the  dynamics  of  a  cantilever  with  a  non-varying  cross  section  are 
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given  by  the  Euler-Bemoulli  beam  equation, 
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where  U (x,  t)  is  the  transverse  beam  defleetion,  n  is  the  mass  per  unit  length,  E  is  Young’s 
modulus,  and  /  is  the  moment  of  inertia  [9].  For  the  ease  of  a  eantilever  where  a  step  foree 
has  been  applied  to  the  tip  at  some  time  in  the  distant  past  the  steady  defleetion  of  the 
eantilever  at  t  =  0  is  given  by 


U{x) 


(6) 


where  L  is  the  length  of  the  eantilever  and  the  appropriate  boundary  eonditions  are  U (0)  = 
17'(0)  =  U”{L)  =  0  and  U'"{L)  =  —Eq/EI.  The  prime  denotes  differentiation  with 
respeet  to  x. 


Similarly,  the  defleetion  of  the  same  eantilever  beam  due  to  the  applieation  of  a  point- 
torque  at  the  eantilever-tip  is  quadratie  in  axial  distance  and  is  given  by 
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where  the  appropriate  boundary  conditions  are  U (0)  =  f/'(0)  =  U"\L)  =  0  and  U'\L)  = 
tq/ EL  The  angle  of  the  cantilever  measured  relative  to  the  horizontal  or  undisplaced 
cantilever  is  then  given  by  tan  0  =  U'{x). 


The  mode  shapes  for  a  cantilevered  beam  are  given  by 


=  —  (cos  kL  +  cosh  k,x)  (cos  kL  —  cosh  kx) 

+  (sin  kL  —  sinh  kx)  (sin  kL  +  sinh  kx)  ,  (8) 

where  n  is  the  mode  number,  and  the  characteristic  frequencies  are  given  by  =  u‘^jj,/EI. 
The  mode  numbers  k  are  solutions  to  1  +  cos  kL  cosh  kL  =  0  [9].  The  initial  cantilever 
displacement  given  by  Eqs.  (6)  and  (7)  can  be  expanded  into  the  beam  modes 

OO 

U{x)  =  ^an^n{x),  (9) 

n=l 

with  mode  coefficients  a„.  The  total  energy  E^  of  the  deflected  beam  is  given  by 

Eb  =  ^  [  U"{xfdx,  (10) 

^  Jo 
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which  is  entirely  composed  of  bending  energy.  The  fraction  of  the  total  bending  energy 
contained  in  an  individual  mode  is  given  by 

(11) 

The  coefficients  for  the  rectangular  cantilever  of  Table  7  are  shown  in  Table  2.  For  the 
case  of  a  force  applied  to  the  cantilever  tip,  97%  of  the  total  bending  energy  is  contained 
in  the  fundamental  mode  and  the  energy  contained  in  the  higher  modes  decays  rapidly 
with  less  than  1%  of  the  energy  contained  in  mode  three.  When  a  point-torque  is  applied 
to  the  same  beam  it  is  clear  that  a  significant  portion  of  the  bending  energy  is  spread  over 
the  higher  modes.  Only  61%  of  the  energy  is  contained  in  the  fundamental  mode  and  the 
decay  in  energy  with  mode  number  is  more  gradual.  The  fifth  mode  for  the  tip-torque  case 
contains  more  energy  than  the  second  mode  for  the  tip-force  case.  Although  we  have  only 
discussed  a  mode  expansion  for  the  rectangular  cantilever,  the  V-shaped  cantilever  will 
exhibit  similar  trends  since  the  transverse  mode  shapes  are  similar  to  that  of  a  rectangular 
beam  [10]. 

The  variation  in  the  energy  distribution  among  the  modes  required  to  describe  the  initial 
deflection  of  the  cantilever  can  be  immediately  connected  to  the  resulting  stochastic  dy¬ 
namics.  For  the  deterministic  calculations  the  initial  displacement  can  be  arbitrarily  set  to 
a  small  value.  In  this  limit  the  modes  of  the  cantilever  beam  are  not  coupled  through  the 
fluid  dynamics.  As  a  result,  the  stochastic  dynamics  of  each  mode  can  be  treated  as  the 
ring  down  of  that  mode  from  the  initial  deflection.  This  indicates  that  the  more  energy  that 
is  distributed  amongst  the  higher  modes  initially  the  more  significant  the  ring  down  and, 
using  the  fluctuation-dissipation  theorem,  the  more  significant  the  stochastic  dynamics. 

The  mode  expansion  clearly  shows  that  the  tip-torque  case  has  more  energy  in  the  higher 
modes.  This  suggests  that  stochastic  measurements  of  the  cantilever  tip-angle  will  have  a 
stronger  signature  from  the  higher  modes  than  measurements  of  cantilever  tip-displacements 
Using  finite  element  simulations  for  the  precise  geometries  of  interest  we  quantitatively 
explore  these  predictions. 

1.2  Computing  the  Driven  Dynamics 

In  order  to  calculate  the  cantilever  dynamics  due  to  an  external  driving  force  we  com¬ 
pute  the  cantilever’s  response  to  an  appropriately  chosen  impulse  in  force.  This  has  been 
done  using  an  impulse  in  velocity  to  explore  the  driven  dynamics  of  cantilevers  beams  of 
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n  bn  (tip-force)  bn  (tip-torque) 


1 

0.97068 

0.61308 

2 

0.02472 

0.18830 

3 

0.00315 

0.06473 

4 

0.00082 

0.03309 

5 

0.00030 

0.02669 

Table  2:  The  fraction  of  the  total  energy  Ei,  contained  in  the  first  five  beam  modes  given  by 
the  coefficients  bn-  The  tip-force  results  are  for  a  rectangular  beam  that  has  been  deflected 
by  the  application  of  a  point  force  to  the  cantilever  tip.  The  tip-torque  results  are  for 
a  rectangular  beam  that  has  been  deflected  by  the  application  of  a  point  torque  to  the 
cantilever  tip.  The  coefficients  clearly  show  that  the  tip-torque  case  has  significantly  more 
energy  contained  in  the  higher  modes. 

(6*f)^'^^(rad) 


(1)  5.6  5.0  X  10-^ 

(2)  20  7.0  X  10-9 


Table  3:  The  magnitude  of  stochastic  fluctuations  in  tip-deflection  and  in  tip-angle  for  the 
rectangular  (1)  and  V-shaped  (2)  cantilevers.  These  values  were  obtained  from  numerical 
simulations  simulations  of  the  beams  in  vacuum. 


varying  geometry,  near  a  solid  wall,  and  including  the  effects  of  higher  modes  of  oscilla¬ 
tion  [11].  In  what  follows  we  focus  upon  the  dynamics  of  the  fundamental  flexural  mode 
and  allow  the  driving  force  to  vary  spatially  given  by. 


Ffibit)  x<i 

0  X  >  ^ 


(12) 


where  x*  =  (L,  5/2,  h/2)  indicates  the  tip  coordinates  where  the  force  is  applied  where 
Fq  is  a  constant  force.  The  time  dependent  displacement  of  the  cantilever  W {x,  t)  due  to 
the  application  of  the  drive  force  is  computed  numerically.  The  power  spectrum  in  terms 
of  cantilever  displacement  is  then  given  by 


Pw{.X,Uj) 


W{x,u) 
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(13) 
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The  power  speetrum  in  terms  of  eantilever  angle  is  found  by  eomputing  the  slope  of 
IV (x,  t)  at  the  region  of  interest  to  yield  0(a;,  t)  and 


Pe{x,uj) 


Q{x,  oj) 
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(14) 


An  advantage  of  this  approaeh  is  that  the  eomplete  speetral  response  over  all  frequency, 
and  for  all  modes,  is  determined  from  a  single  numerical  simulation.  The  alternative  of 
performing  many  simulations  at  different  frequencies  is  computationally  prohibitive  for 
these  systems. 


2  The  role  of  cantilever  geometry 

In  this  section  we  report  on  our  progress  in  quantifying  the  stochastic  dynamics  of  a  can¬ 
tilever  in  fluid  as  a  function  of  the  cantilever  geometry.  The  resulting  fluid-solid  interaction 
problem  is  quite  complex  and  geometric  effects  can  have  a  significant  impact  upon  device 
performance.  We  have  explored  two  important  geometries  that  are  commonly  used:  a  can¬ 
tilever  with  a  rectangular  cross-sectional  area,  and  a  cantilever  with  a  V-shaped  planform 
as  shown  in  Fig.  1.  The  results  of  our  study  have  been  reported  in  detail  in  Refs.  [3,4]. 

2.1  Rectangular  Cantilever 

We  have  performed  deterministic  numerical  simulations  of  the  three-dimensional,  time 
dependent,  fluid-solid  interaction  problem  to  quantify  the  stochastic  dynamics  of  a  rect¬ 
angular  cantilever  immersed  in  water  using  the  thermodynamic  approach  previously  dis¬ 
cussed.  The  deterministic  numerical  simulations  are  done  using  a  finite  element  approach 
that  is  described  elsewhere  [12, 13]. 

The  stochastic  fluctuations  in  cantilever  tip-displacement  for  a  rectangular  cantilever  in 
water  have  been  described  elsewhere  [2, 5, 14, 15].  In  the  following  we  compare  these 
results  with  the  stochastic  dynamics  as  determined  by  the  fluctuations  of  the  cantilever 
tip-angle.  The  geometry  of  the  the  specific  micron  scale  cantilever  we  explore  is  given  in 
Table  7. 

The  autocorrelations  in  equilibrium  fluctuations  follow  immediately  from  the  ring  down 
of  the  cantilever  due  to  the  removal  of  a  step  force  (to  yield  (Mi(O)Mi(t)))  or  step  point- 
torque  (to  yield  (6*i(O)0i(t))).  The  autocorrelations  of  the  rectangular  cantilever  are  shown 
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Figure  1:  Schematics  of  the  two  micron  scale  cantilever  geometries  considered  (not  drawn 
to  scale).  Panel  (a),  A  rectangular  cantilever  with  aspect  ratios  L/h  =  98.5,  w/h  =  14.5, 
and  L/w  =  6.8.  The  cantilever  is  composed  of  silicon  with  density  pc  =  2329  kg/m^ 
and  Youngs  Modulus  E  =  174  GPa.  Panel  (b),  A  V-shaped  cantilever  with  aspect  ratios 
Ljh  =  233,  =  30,  and  L/tu  =  7.8.  The  total  width  between  the  two  arms  normalized 

by  the  width  of  a  single  arm  hh/w  =  10.36.  The  cantilever  planform  is  an  equilateral 
triangle  with  6  =  7r/3.  The  cantilever  is  composed  of  silicon  nitride  with  pc  =  3100kg/m^ 
and  E  =  172GPa.  The  specific  dimensions  for  the  rectangular  and  V-shaped  cantilever 
are  given  in  Table  7. 

in  Fig.  2.  The  magnitude  of  the  noise  is  quantified  by  the  root  mean  squared  tip-angle  and 
deflection  which  is  listed  in  Table  3. 

A  comparison  of  the  autocorrelations  yields  some  interesting  features.  At  short  times 
(6'i(0)6'i(t))  shows  the  presence  of  higher  harmonic  contributions.  This  is  shown  more 
clearly  in  the  inset  of  Fig.  2.  This  further  suggests  that  the  angle  autocorrelations  are  more 
sensitive  to  higher  mode  dynamics. 

The  Fourier  transform  of  the  autocorrelations  yield  the  noise  spectra  shown  in  Fig.  3.  In 
our  notation  the  subscript  of  G  indicates  the  variable  over  which  the  noise  spectrum  is 
measured:  Gg  is  the  noise  spectrum  for  tip-angle  and  Gu  is  the  noise  spectrum  for  tip- 
displacement.  The  equipartition  theorem  of  energy  yields, 

1  k  T 

_  1^  ^  (15) 

T  rG,(i^)dcj  =  (16) 

27r  Jq  kt 

where  k  and  kt  are  the  transverse  and  torsional  spring  constants,  respectively.  The  curves 
in  Fig.  3  are  normalized  using  the  equipartition  result  to  have  a  total  area  of  unity.  Using 
this  normalization  the  area  under  a  peak  is  an  indication  of  the  amount  of  energy  contained 
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Figure  2:  The  normalized  autocorrelation  of  the  rectangular  cantilever  for  tip-deflection 
(solid)  and  tip-angle  (dashed).  (Inset)  A  detailed  view  of  the  autocorrelation  at  short  time 
differences  to  illustrate  the  influence  of  higher  modes  in  the  tip-angle  measurements. 

in  a  particular  mode.  Figure  3  shows  only  the  first  two  modes,  although  the  numerical 
simulations  include  all  of  the  modes  (within  the  numerical  resolution  of  the  finite  element 
simulation).  The  energy  distribution  across  the  first  two  modes  shows  the  significance  of 
the  second  mode  for  the  tip-angle  dynamics. 

Using  a  simple  harmonic  oscillator  approximation  it  is  straight  forward  to  compute  the 
peak  frequency  ujf  and  quality  Q  for  the  cantilever  in  fluid.  Using  a  single  mode  approxi¬ 
mation  yields  the  values  shown  in  Table  4.  As  expected  there  is  a  significant  reduction  in 
the  cantilever  frequency  when  compared  with  the  resonant  frequency  in  vacuum  ojq  and 
the  quality  factor  is  quite  low  because  of  the  strong  fluid  dissipation.  The  values  oiuf  and 
Q  for  tip-angle  and  tip-dispacement  are  nearly  equal.  This  is  expected  since  the  displace¬ 
ments  and  angles  are  very  small,  resulting  in  negligible  coupling  between  the  modes.  Any 
differences  inuf  and  Q  can  be  attributed  to  using  a  single  mode  approximation. 

It  is  useful  to  compare  these  results  with  the  commonly  used  approximation  of  an  oscil¬ 
lating,  infinitely  long  cylinder  with  radius  w/2  [2, 16, 17].  The  cantilever  used  here  has  an 
aspect  ratio  of  L/w  ^  7  and  the  infinite  cylinder  theory  is  quite  good  at  predicting  of  c<;/ 
and  Q. 


9 


Figure  3:  The  noise  spectra  of  stochastic  fluctuations  in  cantilever  tip-angle  (dashed)  and 
tip-deflection  (solid)  for  the  rectangular  cantilever.  The  curves  are  normalized  to  have  the 
same  area,  however  only  the  first  two  modes  are  shown. 

2.2  The  stochastic  dynamics  of  a  V-shaped  cantilever 

We  have  also  explored  the  stochastic  dynamics  of  a  V-shaped  cantilever  in  fluid.  An  inte¬ 
gral  component  of  any  theoretical  model  is  an  analytical  description  of  the  resulting  fluid 
flow  field  caused  by  the  oscillating  cantilever.  The  deterministic  finite  element  simulations 
that  we  performed  yield  a  quantitative  picture  of  the  resulting  fluid  dynamics.  Exploring 
the  flow  fields  further  yields  insight  into  the  dominant  features  that  contribute  to  the  can¬ 
tilever  dynamics. 

As  discussed  earlier,  for  long  and  slender  rectangular  cantilevers  the  flow  field  is  often 
approximated  by  that  of  a  cylinder  of  diameter  w  undergoing  transverse  oscillations.  This 
approach  assumes  that  the  fluid  flow  is  essentially  two-dimensional  in  the  y  —  z  plane 
and  neglects  any  flow  over  the  tip  of  the  cantilever.  Figure  4  (top)  illustrates  this  tip 
flow  for  the  rectangular  cantilever  using  vectors  of  the  fluid  velocity  in  the  x  —  y  plane 
at  z  =  0.  The  figure  is  a  close-up  view  near  the  tip  of  the  cantilever.  It  is  evident  that 
the  flow  over  the  rectangular  cantilever  is  nearly  uniform  in  the  axial  direction  leading 
up  to  the  tip.  However,  near  the  tip  there  is  a  significant  tip  flow  that  decays  rapidly  in 
the  axial  direction  away  from  the  tip.  The  increasing  significance  of  the  tip  flow  as  the 
cantilever  geometry  becomes  shorter  (for  example,  by  simply  decreasing  L)  is  not  certain 
and  remains  an  interesting  open  question.  However,  for  the  geometry  used  here  it  is  clear 
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uif/ujf)  Q 


(1)  0.35  3.34 

(2)  0.36  3.26 

Table  4:  The  peak  frequency  and  quality  factor  of  the  fundamental  mode  of  the  rectangular 
cantilever  determined  by  finite  element  simulations  using  the  thermodynamic  approach. 

(1)  is  computed  using  the  cantilever  tip-displacement  due  to  the  removal  of  a  step  force. 

(2)  is  computed  using  the  cantilever  tip-angle  due  to  the  removal  of  a  point-torque.  The 
frequency  result  is  normalized  by  the  resonant  frequency  in  vacuum  uq.  Using  the  infinite 
cylinder  approximation  with  a  radius  of  w/2  the  analytical  predictions  are  Q  =  3.24  and 
Uf/ujQ  =  0.34. 


that  this  tip-flow  is  negligible  based  upon  the  accuracy  of  the  analytical  predictions  using 
the  two-dimensional  model. 

Figure  4  (bottom)  illustrates  the  tip  flow  for  the  V-shaped  cantilever,  again  by  showing 
velocity  vectors  in  the  x  —  y  plane  dX  z  =  0.  The  shaded  region  indicates  the  part  of 
the  cantilever  where  the  two  arms  have  merged.  To  the  right  of  the  shaded  region  illus¬ 
trates  flow  off  the  tip  and  to  the  left  indicates  flow  that  circulates  back  in  between  the  two 
individual  arms. 

In  order  to  illustrate  the  three-dimensional  nature  of  this  flow,  the  flow  field  in  the  y  —  z 
plane  is  shown  at  two  axial  locations  in  Fig.  5.  Figure  5(top)  is  at  axial  location  x  =  77 ym. 
The  two  shaded  regions  indicate  the  two  arms  of  the  cantilever.  Each  arm  is  generating 
a  flow  with  a  viscous  boundary  layer  (Stokes  layer)  as  expected  from  previous  work  on 
rectangular  cantilevers.  However,  the  Stokes  layers  interact  in  a  complicated  manner  near 
the  center.  It  is  expected  that  as  one  goes  from  the  base  of  the  cantilever  to  the  tip  that 
these  fluid  structures  would  transition  from  non-interacting  to  strongly-interacting. 

Figure  5(bottom)  illustrates  the  flow  field  at  axial  location  x  =  108.8/rm,  the  axial  location 
at  which  the  two  arms  of  the  cantilever  merge  to  form  the  tip  region.  The  length  of  the 
shaded  region  is  therefore  36/im  or  twice  that  of  a  single  arm  shown  in  Fig.  5(top).  For 
this  tip  region  the  flow  field  is  similar  to  what  would  be  expected  of  a  single  rectangular 
cantilever  of  this  width. 

Overall,  it  is  clear  that  the  fluid  flow  field  is  more  complex  for  the  V-shaped  cantilever 
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Figure  4:  The  fluid  flow  near  the  tip  of  the  cantilever  as  illustrated  by  the  velocity  vector 
field  calculated  from  finite  element  numerical  simulations.  A  cross  section  of  the  x  —  y 
plane  at  2;  =  0  is  shown  (see  Fig.  1)  that  is  a  close-up  view  of  the  tip-region.  The  shaded 
region  indicates  the  cantilever  (because  of  the  small  deflections  used  in  the  simulations 
that  cantilever  does  not  appear  to  be  deflected),  (left)  The  flow  field  near  the  tip  of  the 
rectangular  cantilever.  This  flow  field  is  at  t=6ys  and  the  magnitude  of  the  largest  velocity 
vector  shown  is  -0.3  nm/s.  (right)  The  flow  field  near  the  tip  of  the  V-shaped  cantilever. 
This  flow  field  is  at  t=7.2/rs  and  the  magnitude  of  the  largest  velocity  vector  shown  is  -26 
nm/s.  The  shaded  region  indicates  the  tip  region  where  the  two  single  arms  have  merged. 
The  open  region  to  the  left  is  where  the  two  single  arms  have  separated  revealing  the  open 
region  in  the  interior  of  the  V-shaped  cantilever. 

than  for  the  long  and  slender  rectangular  beam.  For  the  V-shaped  cantilever  the  flow  is 
three-dimensional  near  the  tip  region  where  the  two  arms  join  together. 

Central  to  the  flow  field  dynamics  are  the  interactions  of  the  two  Stokes  layers  caused 
by  the  oscillating  cantilever  arms.  The  thickness  of  these  Stokes  layers  are  expected  to 
scale  with  the  frequency  of  oscillation  as  <5^/0  ~  where  a  is  the  half-width  of  the 

cantilever  and  =  uc?  jv  is  a  frequency  based  Reynolds  number  (often  called  the  fre¬ 
quency  parameter).  For  the  relevant  case  of  a  cylinder  of  radius  a  oscillating  at  frequency 
u)  the  solution  to  the  unsteady  Stokes  equations  yields  a  distance  of  approximately  5(5^  to 
capture  99%  of  the  fluid  velocity  in  the  viscous  boundary  layer  [18].  For  a  single  arm  of 
the  V-shaped  cantilever  this  distance  is  nearly  lO/im.  In  comparison,  the  total  distance 
between  the  two  arms  at  the  base  is  125/im.  This  separation  is  large  enough  such  that 
the  two  Stokes  layers  have  negligible  interactions  near  the  base.  However,  as  the  arms 
approach  one  another  with  axial  distance  the  Stokes  layers  overlap  and  eventually  merge 
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Figure  5:  The  fluid  veloeity  veetor  field  at  two  axial  positions  along  the  V-shaped  ean- 
tilever  ealeulated  from  deterministie  finite  element  numerieal  simulations.  Cross  seetions 
of  the  y  —  z  plane  are  shown  (see  Fig.  1),  the  entire  simulation  domain  is  not  shown  and  the 
shaded  region  indieates  the  eantilever.  Both  images  are  taken  at  t=7.2/is  and  the  maximum 
veloeity  veetor  shown  is  -26  nm/s.  (left)  The  y  —  z  plane  at  x  =  77/rm.  The  skewed  width 
of  a  single  arm  of  the  eantilever  in  this  eross-seetion  is  18/im.  The  distanee  separating  the 
two  eantilever  arms  is  36/rm.  (right)  The  y  —  z  plane  at  x  =  108.8/im.  This  is  the  point  at 
whieh  the  two  single  arms  join  to  make  a  eontinuous  eross-seetion  of  width  2w. 

at  the  tip. 

Despite  the  eomplieated  interaetions  of  the  three-dimensional  flow  eaused  by  the  ean¬ 
tilever  tip  and  the  axial  merging  of  the  two  Stokes  layers,  the  V-shaped  eantilever  be¬ 
haves  as  a  damped  simple  harmonie  oseillator.  The  autoeorrelations  in  tip-angle  and  tip- 
displaeement  that  are  found  using  full  finite  element  numerieal  simulations  are  shown  in 
Fig.  6.  It  is  again  elear  that  the  tip-angle  dynamies  have  signifieant  eontributions  from 
the  higher  modes,  see  the  inset  of  Fig.  6.  The  area  normalized  noise  speetra  are  shown  in 
Fig.  7. 

Using  a  simple  harmonie  oseillator  analogy  a  peak  frequeney  and  a  quality  faetor  ean  be 
determined  from  the  first  mode  in  the  noise  speetra  of  Fig.  7.  These  values  are  given  in  the 
first  two  rows  of  Table  6.  The  quality  of  the  eantilever  is  Q  ~  2  and  the  peak  frequeney  is 
redueed  signifieantly,  ujf/ujo  ~  0.2.  eompared  to  the  resonant  frequeney  in  the  absenee  of 
a  surrounding  viseous  fluid. 

It  is  insightful  and  of  praetieal  use  to  determine  the  geometry  of  the  equivalent  reetangular 
beam  that  would  yield  the  preeise  values  of  k,  ujf,  and  Q  ealeulated  for  the  V-shaped 
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Figure  6:  The  normalized  autocorrelation  of  equilibrium  fluctuations  in  the  tip-deflection 
{ui{0)ui{t))  (solid  lined)  and  in  tip-angle  (01(0)6*1  (f))  (dashed-line)  for  the  V-shaped  can¬ 
tilever.  The  inset  shows  a  close-up  of  the  dynamics  for  short  time  differences  to  illustrate 
the  influence  of  the  higher  modes  in  the  tip-angle  measurements. 


cantilever  from  full  finite-element  numerical  simulations.  For  the  rectangular  beam  the 
equations  are  well  known  (c.f.  Ref.  [2])  and  yield  a  unique  value  of  length  L',  width  w', 
and  height  h'  as  shown  below, 


k 

Q 


3EI  Ew'h'^ 


L'3 

7/ 


4L'3 


T"{wkujf) 


(17) 

(18) 


where  the  peak  frequency  is  determined  from  the  maximum  of  the  noise  spectrum, 
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[(1  -  (72(1  +  ToT'{R^u))y  +  {u^nT"{Rou)f] ' 


(19) 


In  the  above  equations  ob  =  cu/cuq  is  the  normalized  frequency,  a  =  0.234  is  a  constant  fac¬ 
tor  to  determine  an  equivalent  lumped  mass  for  a  rectangular  beam,  m/  is  the  equivalent 
mass  of  the  cantilever  plus  the  added  fluid  mass,  7/  is  the  fluid  damping,  F  is  the  hydro- 
dynamic  function  for  an  infinite  cylinder,  F'  is  the  real  part  of  F,  and  F"  is  the  imaginary 
part  of  r.  Equations  (28)-(18)  can  be  solved  to  yield  values  for  the  unknown  geometry  of 
the  equivalent  rectangular  beam  L' ,  w',  and  h'  which  are  given  in  Table  5.  The  equivalent 
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Figure  7:  The  noise  spectra  for  the  V-shaped  cantilever  as  determined  from  the  tip- 
displacement  Gx  (solid  line)  and  from  tip-angle  Ge  (dashed  line).  The  curves  are  nor¬ 
malized  to  have  an  area  of  unity,  with  only  the  first  two  modes  shown. 

beam  is  shorter,  thinner,  and  wider  than  the  V-shaped  cantilever.  Importantly,  the  width  of 
the  equivalent  beam  is  nearly  twice  that  of  a  single  arm  of  the  V-shaped  cantilever. 

L'/L  w'/w  h/h' 


0.8  1.9  0.8 

Table  5 :  The  geometry  of  the  equivalent  rectangular  beam  that  yields  the  exact  values  of  k, 
Uf,  and  Q  for  the  V-shaped  cantilever  that  have  been  determined  from  full  finite-element 
numerical  simulations.  The  length,  width,  and  height  of  the  equivalent  beam  {L',  w',  h') 
are  calculated  using  Eqs.  (28)-(18)  and  are  normalized  by  the  values  of  (L,  w,  h)  for  the 
V-shaped  cantilever  given  in  Table  7. 

These  results  suggest  that  the  parallel  beam  approximation  (PBA)  [19-22]  commonly  used 
to  determine  the  spring  constant  for  a  V-shaped  cantilever  may  also  provide  a  useful  geom¬ 
etry  for  determining  the  dynamics  of  V-shaped  cantilevers  in  fluid.  In  this  approximation 
the  V-shaped  cantilever  is  replaced  by  an  equivalent  rectangular  beam  of  length  L,  width 
2w,  and  height  h  to  yield  a  simple  analytical  expression  for  the  spring  constant.  This  has 
been  shown  to  be  quite  successful  for  V-shaped  cantilevers  that  have  arms  that  are  not  sig¬ 
nificantly  skewed.  The  results  of  using  the  geometry  of  this  approximation  to  determine 
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ojf  and  Q  from  the  two-dimensional  ey Under  approximation  are  shown  on  the  third  row  of 
Table  6.  It  is  elear  that  this  is  quite  aeeurate.  It  is  expeeted  that  these  results  will  remain 
useful  for  eantilever  geometries  that  do  not  deviate  signifieantly  from  that  of  an  equilat¬ 
eral  triangle  as  studied  here.  An  exploration  of  the  breakdown  of  this  approximation  is 
possible  using  the  methods  described  but  is  beyond  the  scope  of  the  current  efforts. 


Ulf/UJQ 

Q 

(1) 

0.21 

1.98 

(2) 

0.22 

2.04 

(L,  2w,  h) 

0.19 

1.98 

Table  6:  The  peak  frequency  and  quality  factor  of  the  fundamental  mode  of  the  V-shaped 
cantilever  determined  by  finite  element  simulations  using  the  thermodynamic  approach. 
(1)  is  computed  using  the  cantilever  tip-displacement  due  to  the  removal  of  a  step  force.  (2) 
is  computed  using  the  cantilever  tip-angle  due  to  the  removal  of  a  point-torque.  The  third 
line  represents  theoretical  predictions  using  the  geometry  of  an  equivalent  rectangular 
beam  given  by  (L,  2w,  h).  The  frequency  result  is  normalized  by  the  resonant  frequency 
in  vacuum  ujq. 


3  The  importance  of  nearby  bounding  surfaces 

In  practice,  the  cantilever  is  never  placed  in  an  unbounded  fluid  and  the  influence  of  nearby 
boundaries  must  be  accounted  for  to  provide  a  complete  description  of  the  dynamics. 
In  many  cases  the  cantilever  is  purposefully  brought  near  a  surface  out  of  experimental 
interest  in  order  to  probe  some  interaction  with  the  cantilever  or  to  probe  the  surface  itself. 
This  can  have  a  significant  impact  upon  the  performance  of  a  device.  Overall,  the  presence 
of  a  boundary  tends  to  increase  the  viscous  dissipation  acting  upon  a  cantilever  causing 
the  quality  factor  of  the  oscillations  to  decrease.  In  our  research  we  have  quantified  this 
precisely  over  a  wide  range  of  experimentally  relevant  conditions.  This  work  has  been 
published  in  Refs.  [3,4]. 

To  specify  our  discussion  we  will  consider  the  situation  depicted  in  Fig.  8  showing  a 
cantilever  a  distance  s  from  a  planar  boundary.  In  the  following  we  study  the  case  where 
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the  cantilever  exhibits  flexural  oscillations  in  the  direction  perpendicular  to  the  boundary. 
However,  we  would  like  to  emphasize  that  our  approach  is  general  and  can  be  used  to 
explore  arbitrary  cantilever  orientations  and  oscillation  directions  if  desired.  The  fluid  is 
assumed  to  be  unbounded  in  all  other  directions.  It  is  well  known  that  the  presence  of  the 
boundary  will  influence  the  dynamics  of  the  cantilever  [23-25].  The  result  is  a  reduction 
in  the  resonant  frequency  and  quality  factor.  This  has  been  described  theoretically  for  the 
case  of  a  long  and  thin  cantilever  of  simple  geometry  where  the  fluid  dynamics  have  been 
assumed  two-dimensional  [26-29]. 


Figure  8:  A  schematic  of  a  cantilever  a  distance  s  away  from  a  solid  planar  surface  (not 
drawn  to  scale).  The  cantilever  undergoes  flexural  oscillations  perpendicular  to  the  sur¬ 
face. 

In  the  following  we  use  the  thermodynamic  approach  with  finite  element  numerical  simu¬ 
lations  to  quantify  the  dynamics  of  the  V-shaped  cantilever  as  a  function  of  its  separation 
from  a  boundary.  We  have  performed  8  simulations  over  a  range  of  separations  from  10  to 
60/im  using  both  the  tip-deflection  and  tip-angle  formulations.  The  noise  spectra  for  these 
simulations  are  shown  in  Fig.  9.  Using  the  insights  from  our  simulations  of  the  V-shaped 
cantilever  in  an  unbounded  fluid  we  expect  the  relevant  length  scale  for  the  fluid  dynamics 
to  be  twice  the  width  of  a  single  arm,  2w.  Using  the  peak  frequency  of  the  V-shaped  can¬ 
tilever  in  unbounded  fluid  yields  a  Stokes  length  5s  =  4.14/rm.  Scaling  the  separation  by 
the  Stokes  length  yields.  2.5  <  s/ 5s  <  15  which  covers  the  range  from  what  is  expected 
to  be  a  strong  influence  of  the  wall  to  a  negligible  influence.  Figure  9  clearly  shows  a 
reduction  in  the  peak  frequency  and  a  broadening  of  the  peak  as  the  cantilever  is  brought 
closer  to  the  boundary.  In  fact,  for  the  smaller  separations  the  peak  is  quite  broad  and  the 
trend  suggests  that  eventually  the  peak  will  become  annihilated  as  the  cantilever  is  brought 
closer  to  the  boundary. 

Using  the  noise  spectra  we  compute  a  peak  frequency  and  a  quality  factor  for  the  funda- 
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Figure  9:  Panel  (a)  The  noise  spectra  Gu  of  stochastic  fluctuations  in  cantilever  tip- 
deflection  for  separations  s  =  10, 12, 15,  20, 40/im.  Panel  (b)  the  noise  spectra  Gg  of 
stochastic  fluctuations  in  cantilever  tip-angle  for  separations  s  =  15,  25,  60/im.  The  spec¬ 
tra  have  been  normalized  by  the  maximum  value  of  Gu  or  Gg.  The  smallest  and  largest 
values  of  separation  are  labeled  with  all  other  values  appearing  sequentially. 

mental  mode  as  a  function  of  separation  from  the  boundary,  which  are  plotted  in  Fig.  10. 
The  horizontal  dashed  line  represents  the  value  of  the  peak  frequency  and  quality  factor  in 
the  absence  of  bounding  surfaces  using  the  two-dimensional  infinite  cylinder  approxima¬ 
tion  [2]  where  the  cylinder  width  has  been  chosen  to  be  2w.  It  is  clear  from  the  results  that 
for  separations  greater  than  s/6s  >  7  the  V-shaped  cantilever  is  not  significantly  affected 
by  the  presence  of  the  boundary.  However,  as  the  separation  decreases  below  this  value 
the  peak  frequency  and  quality  factor  decrease  rapidly. 

The  triangles  in  Fig.  10  represent  the  theoretical  predictions  of  Green  and  Sader  [28, 29] 
using  a  two-dimensional  approximation  for  a  beam  of  uniform  cross-section  that  accounts 
for  the  presence  of  the  boundary.  We  have  used  a  width  of  2w  in  computing  these  the¬ 
oretical  predictions  for  comparison  with  our  numerical  results.  Despite  the  complex  and 
three-dimensional  nature  of  the  flow  field  the  theory  is  able  to  accurately  predict  the  quality 
factor  over  the  range  of  separations  explored.  The  frequency  of  the  peak  for  the  V-shaped 
cantilever  shows  some  deviation  from  these  predictions. 

In  general,  an  increase  in  the  period  of  oscillation  for  a  submerged  object  can  be  attributed 
to  the  mass  of  fluid  entrained  by  the  object  [30].  The  lower  peak  frequency  calculated  for 
the  V-shaped  cantilever  using  a  two-dimensional  solution  indicates  an  over-prediction  of 
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the  mass  loading.  This  can  be  attributed  to  the  three-dimensional  flow  around  the  tip  being 
neglected  for  this  approach.  It  is  reasonable  to  expect  the  cantilever  tip  to  carry  a  smaller 
amount  of  fluid  than  a  section  of  the  beam  body  moving  with  the  same  velocity,  see  Fig.  4. 
The  quality  factor  relates  to  the  ratio  of  the  mass  loading  and  the  viscous  dissipation  and 
is  less  sensitive  to  deviations  incurred  from  the  two-dimensional  approximation.  Despite 
neglecting  three-dimensional  flow  around  the  cantilever  tip,  the  two-dimensional  model 
for  the  fluid  flow  around  the  V-shaped  cantilever  gives  an  accurate  prediction  of  the  peak 
frequency  and  quality  factor. 


Figure  10:  The  variation  of  the  peak  frequency  (panel  (a))  and  quality  (panel  (b))  of  the 
fundamental  mode  of  the  V-shaped  cantilever  in  fluid  as  a  function  of  separation  from 
a  nearby  wall.  Results  calculated  using  tip-deflection  are  circles,  results  using  tip-angle 
are  squares,  and  theoretical  predictions  using  the  results  of  Ref.  [29]  are  triangles.  The 
peak  frequency  and  quality  factor  of  the  fundamental  mode  in  an  unbounded  fluid  are 
cu/cu/  ~  0.19  and  Q  ~  2  and  are  represented  by  the  horizontal  dashed  line.  The  distance 
s  is  normalized  by  the  Stokes  length  5s  where  a  =  w  to  yield  5s  =  4.14/im. 

4  The  dynamics  of  doubly-clamped  beams  in  viscous  fluid : 
tailoring  the  geometry  to  improve  quality  factors 

There  are  numerous  experimental  advantages  to  using  a  doubly-clamped  beam  as  opposed 
to  a  cantilever  in  fluid.  Important  advantages  include  technologies  for  on-chip  sensing  and 
actuation  that  do  not  require  bulky  or  complex  optics.  It  is  also  anticipated  that  the  higher 
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frequency  of  oscillation  of  a  doubly-clamped  beam,  when  compared  to  a  cantilever  of  the 
same  length  and  geometry,  will  lead  to  improved  quality  factors  of  the  oscillations.  In  this 
research  we  have  quantified  this  in  detail  for  a  large  range  of  experimentally  accessible 
parameters.  Our  results  have  been  discussed  in  detail  in  Refs.  [31, 32]. 

4.1  Motivation  for  doubly-clamped  beams 

There  is  a  growing  need  for  fast  and  sensitive  micron  and  nanoscale  sensors  and  actua¬ 
tors  that  operate  in  viscous  fluid  environments.  Many  important  technologies  are  based 
upon  the  dynamics  of  small  elastic  beams  in  fluid  [6,  33-35].  If  an  elastic  beam  is  uni¬ 
formly  reduced  in  size  it  will  become  both  softer  (the  equivalent  spring  constant  is  re¬ 
duced)  and  faster  (the  fundamental  frequency  of  oscillation  increases).  This  advantageous 
trend  is  often  exploited  [35].  However,  in  a  fluid  environment  the  relative  magnitude  of 
viscous  forces  to  inertial  forces  becomes  large  resulting  in  a  dramatic  reduction  in  the 
quality  factor  and  resonant  frequency  of  the  fundamental  mode  of  oscillation.  For  ex¬ 
ample,  the  dynamics  of  a  nanoscale  cantilever  in  water  can  be  overdamped  [1].  Several 
approaches  have  been  proposed  to  overcome  this  difficulty  including  the  use  of  the  higher 
order  beam  modes  [10, 11,33,36-38],  nonlinear  feedback  control  strategies  for  the  ex¬ 
ternal  drive  [39,40],  by  varying  the  cross-sectional  geometry  of  long-thin  cantilevers  that 
are  driven  externally  [11],  and  by  embedding  the  fluid  inside  the  cantilever  while  it  oscil¬ 
lates  in  vacuum  [41].  However,  these  approaches  can  be  difficult  to  implement  in  practice 
and  often  require  sophisticated  measurements  and  control  electronics.  In  addition,  for  the 
strongly  damped  dynamics  under  consideration  here  the  mode  of  actuation  directly  af¬ 
fects  the  resulting  quality  factor  and  resonant  frequency  (c.f.  [2]).  In  many  applications 
a  simpler  tactic  is  desirable  to  overcome  the  strong  viscous  damping.  In  this  paper  we 
explore  the  variation  in  beam  dynamics  as  a  function  of  its  geometry.  In  particular,  we 
quantify  the  stochastic  dynamics  of  doubly-clamped  beams  with  rectangular  cross-section 
for  a  wide  range  of  sizes  and  geometries  including  short  and  wide  beams  that  are  not  well 
described  by  available  analytical  theory.  Using  numerical  simulations  for  the  precise  con¬ 
ditions  of  experiment  we  quantify  the  Brownian  driven  dynamics  of  micron  scale  beams 
in  fluid  and  explore  the  physical  origins  of  the  fluid  dissipation.  These  results  determine 
the  effectiveness  of  tailoring  the  beam  geometry  to  overcome  the  strong  viscous  damping. 

We  calculate  the  stochastic  dynamics  of  the  doubly-clamped  beams,  see  Fig.  11,  using 
the  thermodynamic  approach  discussed  in  detail  in  Refs.  [1,2].  The  approach  requires 
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only  a  single  deterministic  calculation  of  the  fluid  dissipation  that  is  used  to  compute 
the  stochastic  beam  displacement  via  the  fluctuation-dissipation  theorem.  For  a  doubly- 
clamped  beam  the  deterministic  calculation  is  the  ring-down  of  the  beam  due  to  the  re¬ 
moval  of  a  step  point-force  applied  to  the  center  of  the  beam.  We  emphasize  that  the 
only  assumptions  in  this  result  are  that  of  classical  dynamics  and  small  deflections.  Us¬ 
ing  three-dimensional,  time  dependent,  finite  element  simulations  for  the  precise  geome¬ 
tries  of  interest  the  stochastic  dynamics  are  computed.  In  particular,  we  calculate  the 
autocorrelations  and  noise  spectra  of  equilibrium  fluctuations  in  the  beam  displacement. 
The  basic  approach  has  been  validated  both  against  analytics  and  experimental  measure¬ 
ment  [1,2,26,27,42]  and  also  used  to  study  the  fluid-coupled  motion  of  two  atomic  force 
microscope  cantilevers  [43]  and  two  nanoscale  cantilevers  [1,2].  In  many  situations  of 


Figure  11:  A  schematic  of  a  doubly-clamped  beam  used  in  the  numerical  simulations 
with  length  L,  width  w,  and  height  h  with  uniform  rectangular  cross-section,  (a)  The 
X  —  z  plane  of  the  beam  in  fluid.  The  beam  is  supported  by  a  rigid  support  of  width 
w  on  each  side  for  the  short-wide  geometries  to  minimize  the  effects  of  the  bounding 
side  walls.  The  beam  is  light  grey  and  the  two  rigid  supports  are  darker  grey,  (b)  The 
y  —  z  plane  of  the  beam  illustrating  the  rectangular  cross-section.  In  our  simulations  the 
beam  is  immersed  in  room  temperature  water  and  we  compute  the  stochastic  dynamics  of 
the  fundamental  flexural  mode  driven  by  Brownian  motion.  In  the  following  figures  the 
vertical  displacement  of  the  beam  at  a;  =  L/2  is  referred  to  as  Zi{t)  for  thermally  induced 
fluctuations  and  Zi{t)  for  the  deterministic  ring  down  simulations. 

technological  and  scientific  interest,  such  as  atomic  force  microscopy,  the  elastic  beams 
are  long  and  thin  L  3>  te  3>  where  L  is  the  length,  w  is  the  width,  and  h  is  height  of 
the  beam.  The  fluid-solid  interaction  problem  describing  the  motion  of  a  waving  beam  in 
fluid  is  very  difficult  with  analytical  solutions  available  only  under  idealized  conditions 
such  as  simple  beam  geometries  and  for  small  deflections  [16,  17,44,45].  In  the  limit 
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Figure  12:  (a)  The  predicted  variation  of  the  quality  factor  Q  for  the  stochastic  displace¬ 
ment  of  a  beam  immersed  in  a  viscous  fluid  with  respect  to  the  nondimensional  frequency 
parameter  Rq  and  mass  loading  parameter  Tq.  (b)  The  predicted  variation  of  the  resonant 
frequency  in  fluid  ojf  with  respect  to  Rq  and  Tq.  In  both  panels  five  curves  are  shown 
for  Tq  =  0.5, 1,  2, 4,  8.  The  bounding  two  curves  are  labeled  with  the  remaining  curves 
in  sequential  order.  Rq  is  evaluated  at  the  resonant  frequency  of  the  beam  in  vacuum. 
The  quality  Q  is  determined  by  evaluating  Eq.  (26)  at  Uf  where  cof  is  the  frequency  that 
maximizes  Eq.  (20). 

of  small  beam  displacements,  a  two-dimensional  approximation  for  the  fluid  flow  over 
the  beam  is  often  used  to  determine  the  force  interactions  with  an  Euler-Bernoulli  beam. 
This  approach  has  been  very  successful  in  predicting  the  resulting  beam  dynamics  in  a 
viscous  fluid  [16].  Eurthermore,  it  has  been  shown  that  replacing  the  rectangular  beam 
cross-section  with  that  of  a  cylinder  of  diameter  equal  to  the  width  w  yields  small  errors 
on  the  order  of  several  percent  [17].  The  flow  field  generated  by  an  oscillating  cylinder 
is  well  known  as  well  as  the  forces  acting  on  the  surface  of  the  cylinder  [30,46].  These 
approximations  have  led  to  insightful  analytical  expressions  describing  the  stochastic  dy¬ 
namics  of  beams  in  fluid  [1, 16].  However,  the  validity  and  accuracy  of  these  expressions 
remain  unclear  for  the  finite  beam  geometries  often  used  in  experiment. 

In  the  limit  of  a  long  and  thin  beam,  small  displacements,  and  using  the  two-dimensional 
approximation  of  an  oscillating  cylinder  for  the  fluid  flow,  the  noise  spectrum  of  equilib¬ 
rium  fluctuations  in  displacement  of  the  beam  measured  at  a;  =  L/2  for  the  fundamental 
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mode  is  given  by  [1], 


^ksT  1 _ TouVi{RQU}) _ 

k  uo  [(1  -  cD2  (1  +  ToTr{RoCo))f  +  {u^ToT,{Rou)f]  ’ 


(20) 


where  u  is  the  frequeney  of  oseillation,  u  =  cu/c^o  is  the  redueed  frequency,  uq  is  the 
resonant  frequency  of  the  fundamental  mode  in  vacuum,  Rq  is  the  frequency  parameter 
evaluated  at  uo,  r(c<;)  is  the  hydrodynamic  function,  Tq  is  the  mass  loading  parameter, 
ks  is  Boltzmann’s  constant,  T  is  the  temperature,  and  k  is  the  spring  constant  for  the 
fundamental  mode.  The  frequency  parameter  is. 
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and  is  a  frequency  based  Reynolds  number  representing  the  ratio  of  local  inertia  to  viscous 
forces  where  v  is  the  kinematic  viscosity  of  the  fluid.  In  our  notation,  the  frequency 
parameter  R  is  evaluated  at  arbitrary  frequency  uj,  and  Rf  is  evaluated  at  cu/.  The  mass 
loading  parameter  is 


Tn  = 
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4  pbh 

and  represents  the  ratio  of  the  mass  of  a  cylinder  of  fluid  with  radius  w/2to  the  actual  mass 
of  the  beam  where  p/  is  the  density  of  the  fluid,  and  pb  is  the  density  of  the  beam.  The 
hydrodynamic  function  for  an  oscillating  cylinder  in  a  viscous  fluid  is  given  by  [30,46], 


r(w) 
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where  Ki  and  Kq  are  Bessel  functions,  T^  and  Tj  are  the  real  and  imaginary  parts  of  T, 
respectively  and  i  = 


The  dynamics  of  a  beam  in  fluid  are  not  precisely  equivalent  to  that  of  a  damped  simple 
harmonic  oscillator.  For  example,  both  the  mass  and  damping  are  frequency  dependent. 
The  mass  of  the  entrained  fluid  plus  the  mass  of  the  beam  is 


mf{u)  =  me  (1  +  Tor^(i?oa;))  (24) 

where  rrie  =  arub  is  the  equivalent  mass  of  the  beam  such  that  the  kinetic  energy  of  this 
mass  is  equal  to  that  of  the  fundamental  mode  and  rub  =  pbLwh  is  the  mass  of  the  beam. 
For  the  fundamental  flexural  mode  of  a  doubly-clamped  beam  a  =  0.396.  The  viscous 
damping  is 

7/(ct;)  =  rricyi^e^V  i{RQU))  (25) 
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where  mcyi,e  =  arricyi  is  the  equivalent  mass  of  a  eylinder  of  fluid  with  diameter  equal 
to  te.  As  the  frequency  of  oscillation  increases  the  magnitude  of  m/  decreases  and  the 
magnitude  of  7/  increases. 

The  simple  harmonic  oscillator  approximation  is  convenient  to  define  commonly  used 
diagnostics  such  as  the  quality  factor  Q  and  resonant  frequency  of  the  beam  in  fluid  cu/.  As 
a  result  of  the  frequency  dependent  mass  and  damping,  the  fundamental  peak  of  the  noise 
spectra  is  not  well  approximated  as  a  Lorentzian  for  these  strongly  damped  oscillators  and 
care  must  be  taken  when  determining  Q  and  Uf.  The  resonant  frequency  in  fluid  u /  will  be 
defined  to  be  the  frequency  which  maximizes  the  noise  spectrum  in  Eq.  (20).  The  quality 
factor  Q  is  then  defined  as  the  ratio  of  energy  stored  by  the  potential  and  kinetic  energy  of 
the  beam  and  fluid  to  the  energy  dissipated  by  viscosity  per  oscillation  when  evaluated  at 
00 f.  This  yields 

r,  _  ^/(^/)^/  _  +  TrjRoUf) 

vi^f)  URoCof)  •  ^  ^ 

Given  values  of  the  nondimensional  parameters  Rq  and  Tq,  Eqs.  (20)  and  (26)  directly 
yield  the  analytical  predictions  for  ojj  and  Q.  The  variation  of  Q  and  ujf  with  Rq  and 
To  are  shown  in  Eig.  12  over  a  large  range  of  parameters.  The  quality  factor  increases 
significantly  as  the  frequency  of  oscillation  is  increased  and  also  increases  as  the  mass 
loading  decreases.  The  resonant  frequency  of  the  beam  when  placed  in  fluid,  uJflujQ,  also 
increases  with  frequency  of  oscillation  and  with  a  reduction  in  mass  loading.  The  increase 
of  c<;//c<;o  with  respect  to  Rq  is  very  rapid  for  Rq  <  20  with  only  small  changes  for  higher 
frequencies,  while  the  dependence  upon  Tq  results  in  a  nearly  uniform  increase  over  the 
range  shown.  It  is  typical  for  i?o  ~  1  and  Tq  ~  1  for  many  proposed  microscale  applica¬ 
tions  in  water.  In  this  case  the  analytics  predict  strongly  damped  dynamics  with  Q  ~  2. 
Eor  applications  that  require  a  distinct  peak  to  be  measured  this  presents  a  significant 
challenge. 

Using  Euler-Bernoulli  beam  theory  [9]  for  a  doubly-clamped  beam  these  expressions  can 
be  written  as  a  function  of  geometry  {L,w,h)  which  are  often  the  experimentally  relevant 
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parameters  rather  than  Rq  and  Tq.  The  relevant  expressions  are, 


(27) 

(28) 
(29) 


where  E  is  the  Youngs  modulus.  These  expressions  for  Tq  and  Rq  together  with  Fig.  12 
suggest  that  Q  and  uj  inerease  by  reducing  the  length,  increasing  the  width,  or  increasing 
the  height  of  the  beam.  However,  the  precise  improvement  is  not  clear  since  the  avail¬ 
able  theoretical  predictions  are  only  for  long  and  slender  beams.  In  light  of  this  we  have 
performed  full  time-dependent  and  three-dimensional  finite  element  numerical  simula¬ 
tions  [13]  of  a  wide  range  of  geometries  to  determine  precisely  the  stochastic  dynamics. 


To  compute  the  stochastic  dynamics  of  the  beams  we  use  the  approach  discussed  in 
Refs.  [1,2]  and  provide  only  the  essential  details  necessary  for  our  discussion.  The  auto¬ 
correlation  of  equilibrium  fluctuations  in  beam  displacement  are  given  by  the  deterministic 
ring-down  of  the  beam  to  the  removal  of  a  point  step  force  applied  at  a;  =  L/2  given  by 


m 


Fq  for  f  <  0 
0  for  t  >  0 


(30) 


where  t  is  time,  and  Fq  is  the  magnitude  of  the  force.  The  value  of  Fq  is  chosen  for  each 
simulation  such  that  the  beam  deflections  remain  small  and,  in  this  case,  the  results  are 
independent  of  its  specific  value.  The  autocorrelation  of  equilibrium  fluctuations  in  beam 
displacement  is  then  given  by, 

{z,{0)z,{t))  =  ksT^.  (31) 

-To 

We  use  lower  case  zi  to  indicate  stochastic  displacement,  and  upper  case  Zi  to  indicate  the 
deterministic  ring-down  measured  at  the  center  of  the  beam  x  =  Lj^.  The  noise  spectrum 
of  fluctuations  in  beam  displacement  is  given  by. 


G{(jj)  =  A  /  {zi{t)zi{R))  cos  {ojt)  dt. 


(32) 


Jo 

The  noise  spectrum  is  used  to  determine  Uf  and  Q  for  the  numerical  results.  The  resonant 
frequency  cuf  is  the  frequency  maximizing  G{u)  and  the  quality  is  given  by 
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The  right  hand  side  of  Eq.  (33)  is  found  using  mf{ujf)  =  ^/kJuJf  and  using  the  peak  value 
of  the  noise  speetrum  G{ujf  )  to  determine  the  damping.  The  error  in  using  the  bulk  mode 
spring  eonstant,  as  opposed  to  the  dynamie  spring  eonstant  for  the  fundamental  mode,  is 
small  and  on  the  order  of  several  pereent. 

In  summary,  the  numerieal  proeedure  is  the  following:  (i)  Compute  Z{t)  from  a  deter- 
ministie  simulation  of  the  ring  down  of  the  beam  due  to  the  removal  of  a  step  foree;  (ii) 
Compute  the  autoeorrelation  of  equilibrium  fluetuations  in  displaeement  using  Eq.  (31); 
(iii)  Compute  the  noise  speetrum  using  Eq.  (32);  (iv)  Caleulate  diagnosties:  uj  is  the 
frequeney  that  maximizes  the  noise  speetrum,  and  Q  is  found  from  Eq.  (33). 

We  have  performed  extensive  numerieal  tests  on  doubly-elamped  beams  in  vaeuum  and 
in  fluid  to  ensure  the  aeeuraey  of  our  ealeulations  [31].  Eor  eaeh  geometry  explored  we 
have  eondueted  numerieal  simulations  over  a  range  of  spatial  and  temporal  diseretiza- 
tions  to  ensure  the  eonvergenee  of  our  reported  values  for  the  quality  faetor  and  resonant 
frequeney  of  the  fundamental  mode  in  fluid.  The  required  spatial  resolution  depended  sig- 
nifieantly  upon  the  geometry  explored  with  the  short  and  wide  beam  geometries  requiring 
higher  spatial  resolution.  Typieally,  we  found  that  a  time  step  At  <  P/15  was  suffieient 
where  P  is  the  period  of  the  fundamental  mode  in  vaeuum.  We  have  also  been  eareful  to 
ehoose  the  size  of  the  overall  simulation  domain  to  be  large  enough  sueh  that  the  bounding 
walls  do  not  affeet  the  results.  In  our  results,  the  bounding  walls  are  always  a  distanee  of 
15(is  or  greater  from  the  beam  surfaee  where  5s  =  {u/ufY^'^  is  the  Stokes  length  for  the 
fundamental  mode  in  fluid. 

As  the  baseline  geometry  we  eonsider  a  doubly-elamped  beam  with  length  L'  =  15/im, 
width  w'  =  0.4/im,  and  height  h'  =  O.l/rm.  This  geometry  is  similar  to  what  has  been  re- 
eently  used  in  experiments  demonstrating  thermoelastie  aetuation  in  vaeuum  and  air  [33]. 
Here,  we  are  interested  in  the  beam  dynamies  in  a  viseous  fluid  and  use  water.  This  geom¬ 
etry  is  referred  to  as  ease  (1)  in  Table  7  and  we  eonsider  seven  additional  geometries  whieh 
are  ehosen  as  systematie  variations  of  the  baseline  geometry  (L',  w\  h').  Also  shown  in 
Table  7  are  the  aspeet  ratios  for  the  different  geometries  to  give  an  idea  of  the  range  of 
geometries  used  and  also  to  give  some  indieation  of  the  deviation  from  the  ideal  ease  of  a 
long  and  thin  beam  used  in  analytieal  predietions. 

Table  8  illustrates  the  deviations  in  geometry  when  eompared  with  the  baseline  geometry 
of  ease  1 .  Also  ineluded  are  the  beam  properties  that  ean  be  determined  independent  of 
the  fluid  dynamies  whieh  inelude  the  bulk  spring  eonstant  k,  the  frequeney  parameter  in 


26 


vacuum  Rq,  and  the  mass  loading  parameter  Tq.  We  have  used  finite  element  numerieal 
simulations  of  the  beams  in  vaeuum  to  determine  the  numerieal  values  of  k  and  uq  for 
all  of  the  geometries  eonsidered.  Given  this  information  one  ean  use  the  analytieal  ex¬ 
pressions  to  prediet  Q  and  cc/  whieh  is  illustrated  in  Fig.  12.  From  Table  8  it  is  elear 
that  over  four  orders  of  magnitude  of  spring  eonstant,  over  three  orders  of  magnitude  of 
frequeney  parameter,  and  over  a  single  order  of  magnitude  of  the  mass  loading  parameter 
are  eonsidered  by  the  ehosen  variations  in  geometry. 


Case 

L  [yum] 

w  [/im] 

h  [/im] 

L/w 

L/h 

w/h 

(1) 

15 

0.4 

0.1 

37.5 

150 

4 

(2) 

15 

0.8 

0.1 

18.75 

150 

8 

(3) 

15 

1.2 

0.1 

12.5 

150 

12 

(4) 

15 

0.4 

0.2 

37.5 

75 

2 

(5) 

15 

0.4 

0.3 

37.5 

50 

1.33 

(6) 

5 

0.4 

0.1 

12.5 

50 

4 

(7) 

1 

0.4 

0.1 

2.5 

10 

4 

(8) 

0.4 

0.4 

0.1 

1 

4 

4 

Table  7:  The  eight  geometries  of  doubly-elamped  beams  used  in  the  numerieal  simula¬ 
tions.  Case  (1)  is  the  baseline  geometry  and  the  remaining  eases  are  variations  of  this 
geometry.  The  beam  aspeet  ratios  are  L/w,  L/h,  and  w/h.  The  horizontal  lines  sepa¬ 
rate  the  different  studies  performed:  the  baseline  geometry,  variation  in  width,  variation 
in  height,  and  variation  in  length.  The  beams  are  eomposed  of  silieon  with  Young’s  mod¬ 
ulus  E  =  210  GPa,  density  pb  =  3100  kg/m^  and  the  fluid  is  water  with  pf  =  997 
kg/m^,  7]  =  8.56  X  10“^  kg/ms.  All  simulations  are  performed  at  room  temperature  with 
T  =  300K. 

We  first  quantify  the  stoehastie  dynamies  of  the  baseline  geometry.  The  numerieal  results 
for  the  autoeorrelation  of  equilibrium  fluetuations  in  beam  displaeement  are  shown  in 
Fig.  13(a),  and  the  noise  speetrum  is  shown  in  Fig.  13(b).  In  eaeh  figure  the  baseline 
geometry  is  labeled  w'. 

The  autoeorrelation  eurves  are  normalized  using  k/ksT  where  the  value  of  k  for  eaeh 
ease  is  given  in  Table  8.  The  noise  speetra  have  been  normalized  using  the  peak  value 
G{(jjf  ).  These  figures  illustrate  that  the  dynamies  of  this  mieron  seale  beam  in  water 
are  strongly  damped.  The  value  of  the  quality  and  resonant  frequeney  in  fluid  using  our 
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Case 

L/L' 

w/w' 

h/h! 

k  [N/m] 

Rq 

To 

(1) 

1 

1 

1 

0.40 

1.08 

1.01 

(2) 

1 

2 

1 

0.83 

4.55 

2.02 

(3) 

1 

3 

1 

1.21 

10.19 

3.03 

(4) 

1 

1 

2 

3.19 

2.11 

0.51 

(5) 

1 

1 

3 

10.82 

3.19 

0.34 

(6) 

1/3 

1 

1 

10.92 

10.20 

1.01 

(7) 

1/15 

1 

1 

1178.95 

231.55 

1.01 

(8) 

1/37.5 

1 

1 

11413.04 

1146.08 

1.01 

Table  8:  The  geometry  variations  with  respect  to  the  baseline  geometry  given  by  case  (1) 
with  (L',  w\  h').  Cases  (2)  and  (3)  explore  increasing  width,  cases  (4)  and  (5)  explore  in¬ 
creasing  thickness,  and  cases  (6)-(8)  explore  decreasing  length.  Also  shown  are  the  spring 
constant  k,  the  frequency  based  Reynolds  number  in  vacuum  Rq,  and  the  mass  loading 
parameter  Tq.  The  values  of  k  and  Rq  are  determined  using  finite  element  simulations. 

numerical  results  are  given  in  Table  9  and  are  Q  =  0.80  and  uJflujQ  =  0.22.  Also  shown 
are  the  predictions  from  analytics  using  Eqs.  (20)  and  (26)  which  yield  =  0.68  and 
oOf  / 00 Q  =  0.22.  The  analytical  predictions  are  quite  accurate  for  the  frequency  drop  while 
under  predicting  the  quality  factor  for  this  geometry. 

Next  we  consider  the  variation  in  the  stochastic  dynamics  of  the  beam  as  a  function  of  the 
beam  width.  In  particular,  we  double  and  triple  the  beam  width  w  while  holding  L  and 
h  constant.  For  increasing  width  the  frequency  parameter  increases  as  Rq  ~  w'^  while 
the  mass  loading  parameter  increases  as  Tq  ~  w.  This  has  the  effect  of  increasing  the 
fluid  inertia  while  simultaneously  increasing  the  mass  loading.  These  two  counteracting 
effects  suggest  the  increase  in  Q  and  oOf  will  only  be  moderate.  The  autocorrelations 
and  noise  spectra  from  numerical  simulations  are  shown  in  Fig.  13.  The  autocorrelation 
results  exhibit  both  positively  and  negatively  correlated  results  as  expected  with  the  dy¬ 
namics  becoming  more  underdamped  as  the  width  is  increased.  The  noise  spectra  clearly 
illustrate  that  the  peak  value  shifts  to  higher  frequency  and  that  the  peak  itself  becomes 
sharper  as  the  width  increases.  For  case  1,  the  noise  spectra  has  significant  contributions 
at  low  frequency  whereas  for  case  3  the  noise  spectra  has  become  more  symmetric  with  a 
Forentzian  shape. 

Values  for  the  quality  and  resonant  frequency  in  fluid  from  our  numerical  results  are  given 
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Case 

R 

Ufjuo 

Q 

Q/Q' 

Uf/ul 

(1) 

0.23 

0.22 

0.80 

1.0 

0.22 

0.68 

(2) 

1.08 

0.24 

1.03 

1.29 

0.27 

1.05 

(3) 

2.91 

0.29 

1.36 

1.71 

0.28 

1.31 

(4) 

1.00 

0.46 

1.28 

1.61 

0.50 

1.35 

(5) 

2.11 

0.66 

2.01 

2.52 

0.66 

2.13 

(6) 

4.80 

0.47 

1.57 

1.97 

0.51 

2.04 

(7) 

164.46 

0.71 

6.13 

7.69 

0.67 

9.22 

(8) 

885.19 

0.77 

5.90 

7.40 

0.69 

20.24 

Table  9:  The  stochastic  dynamics  of  the  beams  in  fluid.  Shown  is  the  frequency  based 
Reynolds  number  in  fluid  R,  the  reduction  in  the  resonant  frequency  Ufjuo,  and  the  quality 
factor  Q.  Also  shown  is  the  improvement  of  the  quality  with  respect  to  that  of  case  (1) 
given  by  Q'  =  0.8.  and  uifjujQ  are  the  results  predicted  from  analytical  theory  using 
Eqs.  (20)  and  (26). 

in  Table  9.  When  compared  to  the  quality  for  the  baseline  geometry  Q',  the  increase  in 
quality  is  Q/Q'  =  1.29  for  doubling  the  width,  and  Q/Q'  =  1.71  for  tripling  the  width. 
The  quality  increases  with  increasing  width  however  the  magnitude  of  the  quality  is  small 
indicating  that  the  beam  dynamics  remain  strongly  damped.  The  increase  in  the  value  of 
ojf/ojQis  slightly  less  than  what  is  found  for  Q.  A  comparison  of  our  numerical  values  of 
ojf  and  Q  with  the  analytical  predictions  of  Eqs.  (20)  and  (26)  are  shown  in  Eig.  14.  The 
circles  are  the  results  from  our  numerical  simulations  and  the  dashed  line  is  the  analytical 
prediction  where  ^  =  w/w'  and  w'  is  the  width  of  the  baseline  geometry.  It  is  clear  that 
that  the  analytical  predictions  remain  quite  accurate  over  this  range.  This  includes  case  3 
where  L/w  ~  w/h  ~  12  and  L/h  3>  1.  Eigure  14  indicates  that  the  magnitude  of  the 
increase  in  quality  with  increasing  width  is  quite  moderate.  Eurthermore,  the  increase  in 
cu/  is  quite  small  and  becomes  nearly  flat  at  cuy/cuo  ~  0.24  for  >  2. 

Next  we  consider  the  variation  in  beam  dynamics  as  the  height  is  increased.  We  consider 
the  cases  where  h  is  doubled  and  tripled  while  the  L  and  w  are  held  constant.  As  the 
height  is  increased  the  frequency  parameter  increases  as  i?o  ~  ^  whereas  the  mass  loading 
parameter  decreases  as  Tq  ~  h~^ .  These  two  effects  both  contribute  to  increasing  Q  and 
ojf.  The  normalized  autocorrelations  and  noise  spectra  from  our  numerical  simulations 
are  shown  in  Eig.  15.  The  results  clearly  indicate  an  increasing  value  of  both  ojf  and  Q 
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Figure  13:  The  autocorrelations  and  noise  spectra  of  equilibrium  fluctuations  in  beam 
displacement  as  a  function  of  beam  width  from  numerical  results:  case  1  (w'),  case  2 
(2w'),  case  3  (3w').  (a)  The  autocorrelations,  the  results  have  been  normalized  using 
k/{kBT)  for  each  case,  (b)  The  noise  spectra,  the  results  have  been  normalized  by  the 
peak  value  for  each  case,  G{ujf). 

and  numerical  values  are  given  in  Table  9.  The  relative  increase  in  quality  isQ/Q'  =  1.61 
when  the  height  is  doubled,  and  Q/Q'  =  2.52  when  the  height  is  tripled.  The  increase  in 
ujf  I Uq  follows  a  similar  trend. 

A  comparison  of  our  numerical  results  with  the  predictions  of  theory  is  shown  in  Fig.  14 
using  ^  =  h/h' .  The  square  symbols  are  the  numerical  results  and  the  solid  line  is  the  an¬ 
alytical  prediction.  It  is  clear  that  the  increases  in  Q  and  cu/  are  much  larger  for  variations 
in  height  when  compared  to  what  was  found  for  increases  in  beam  width.  The  analytical 
predictions  remain  quite  accurate  and  insightful  over  the  range  of  aspect  ratios  explored 
by  varying  the  beam  height.  We  highlight  that  this  includes  case  5  where  ~  1. 

The  last  case  we  consider  is  decreasing  the  beam  length  while  holding  the  width  and 
height  constant.  In  this  case  the  frequency  parameter  increases  rapidly  as  Rq  ~ 
whereas  To  remains  constant.  The  autocorrelations  and  noise  spectra  are  shown  in  Fig.  16 
which  illustrate  a  significant  increase  in  resonant  frequency  and  quality.  From  Fig.  16(a) 
the  results  for  the  most  extreme  geometry  explored,  T'/37.5,  clearly  show  the  influence 
of  higher  harmonics.  The  numerical  values  of  Q  and  ojf  from  our  numerical  results  are 
given  Table  9.  For  case  8  where  LjL'  =  37.5  the  increase  in  quality  is  QjQ'  =  7.4 
and  the  reduction  of  the  resonant  frequency  when  compared  to  its  value  in  vacuum  is 
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Figure  14:  (color)  Comparison  of  numerical  results  with  analytical  predictions  for  Q  and 
a;/  as  a  function  of  width  w  and  height  h.  The  circles  (blue)  are  for  increasing  width  and 
the  squares  (red)  are  for  increasing  height.  The  solid  line  is  the  analytical  prediction  for 
increasing  width  and  the  dashed  line  is  the  analytical  prediction  for  increasing  height.  To 
place  all  values  on  a  single  plot  ^  =  w/w'  for  the  varying  width  results  and  ^  =  h/h'  for 
the  varying  height  results. 

LVfltvo  =  0.77  indicating  significant  changes  are  possible  by  changing  the  beam  length. 

The  analytical  predictions  given  in  Table  9  show  significant  deviations  from  our  numerical 
results.  This  is  also  illustrated  in  Fig.  17  where  the  triangles  are  the  numerical  results  and 
the  solid  lines  are  the  analytical  predictions.  For  case  6  (L/w  =  12.5)  the  analytical 
predictions  are  quite  accurate.  However,  for  case  7  {L/w  =  2.5),  and  case  8  {L/w  =  1) 
the  analytical  predictions  over  predict  Q  and  under  predict  Uf/uQ.  The  approximation 
of  using  the  fluid  flow  from  an  infinite  two-dimensional  oscillating  cylinder  is  no  longer 
well  justified.  The  numerical  results  suggest  the  presence  of  additional  modes  of  fluid 
dissipation  that  are  not  captured  in  the  two-dimensional  theory. 

To  explore  this  further  we  quantify  the  fluid  motion  around  the  beam  in  the  determinis¬ 
tic  numerical  simulations  where  the  beam  rings  down  upon  the  removal  of  a  step  force. 
Figures  18  and  19  illustrate  the  magnitude  of  the  fluid  velocity  in  the  transverse  and 
axial  Ux  directions,  respectively.  The  velocities  are  plotted  along  a  line  beginning  at 
(0.0,  0.0,  0.01/xm)  and  ending  at  (L,  0.0,  O.Ol/rm)  for  all  cases.  In  our  notation  the  fluid 
velocity  in  the  (x,  y,  z)  directions  is  {Ux,  Uy,  Uz),  see  Fig.  1 1  for  the  definition  of  the  coor¬ 
dinate  directions  {x,y,  z).  The  velocities  are  shown  at  the  time  when  the  velocity  of  the 
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Figure  15:  The  autocorrelations  and  noise  spectra  as  a  function  of  beam  height  from  nu¬ 
merical  results:  case  1  (h'),  case  4  (2h'),  case  5  (3h').  (a)  The  normalized  autocorrelations, 
(b)  The  normalized  noise  spectra. 

beam  is  at  its  maximum  value  which  occurs  when  the  center  of  the  beam  crosses  z  =  0 
the  first  time  during  its  ring  down.  The  maximum  value  of  at  this  time  is  labeled  Mq 
and  is  used  to  normalize  both  the  transverse  and  axial  velocities.  The  axial  direction  is 
normalized  by  the  length  so  that  all  cases  can  be  represented  on  the  same  figure. 

Figure  18  shows  the  transverse  fluid  velocity  for  cases  1-8.  The  baseline  geometry 
(dashed  line)  and  cases  2-7  (solid  lines)  collapse  onto  a  single  curve  with  a  shape  similar 
to  that  of  the  fundamental  mode  of  a  doubly-clamped  beam.  Case  8  differs  significantly 
with  a  much  sharper  peak  indicating  that  its  dynamics  are  quite  different  which  is  expected 
since  this  geometry  is  substantially  different  than  the  others. 

Figure  19(a)  illustrates  the  normalized  axial  velocities  as  the  beam  width  and  height 
are  varied.  In  the  approximation  of  a  two-dimensional  flow  the  axial  velocity  is  identically 
zero  and  any  deviations  from  this  in  the  numerical  results  indicate  fluid  dynamics  not 
considered  in  the  analytical  predictions.  The  bimodal  shape  of  the  curves  is  expected 
from  the  symmetry  of  the  fundamental  mode.  For  x  >  L/2  the  axial  fluid  velocity  is 
positive  and  for  a;  <  L/2  it  is  negative.  The  baseline  geometry  is  shown  as  the  dashed 
line  and  has  a  negligible  axial  fluid  velocity.  At  its  maximum  value  it  is  only  ~  2.5% 
of  the  maximum  transverse  velocity  uq.  A  similar  trend  is  found  for  cases  2-5.  As  the 
width  or  height  is  increased  the  relative  magnitude  of  the  axial  velocities  increase.  It  is 
expected  that  if  larger  values  of  the  width  or  height  were  computed  the  axial  velocities 
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Figure  16:  The  autocorrelations  and  noise  spectra  as  a  function  of  beam  length  from  nu¬ 
merical  results:  case  1  {L'),  case  6  (T'/3),  case  7  (T'/5),  case  8  (L'/ST.S).  (a)  The  nor¬ 
malized  autocorrelations,  (b)  The  normalized  noise  spectra. 

would  become  significant  and  at  this  point  the  analytical  predictions  would  show  large 
deviations. 

Figure  19(b)  shows  the  relative  value  of  the  axial  velocities  as  the  length  of  the  beam 
is  decreased.  The  baseline  geometry  is  included  for  reference  as  the  dashed  line.  It  is 
clear  that  the  axial  velocities  are  now  quite  significant  and  range  from  10%  to  40%  of 
Wo-  The  axial  velocities  do  not  vanish  at  a;  =  0,  L  for  cases  7  and  8  because  these  beams 
are  held  by  rigid  supports,  see  Fig.  11,  and  the  lateral  side  walls  of  the  numerical  domain 
are  distant.  The  axial  velocities  result  in  fluid  dissipation  not  accounted  for  in  the  two- 
dimensional  theory  and  contribute  significantly  to  the  lower  values  of  Q  found  in  the 
numerical  simulations.  Furthermore,  Uf  from  the  numerical  simulations  are  larger  than  the 
analytical  predictions.  The  added  mass  in  the  simulations  are  smaller  than  the  predicted 
values  and  this  reduction  is  a  direct  result  of  the  three-dimensionality  of  the  fluid  flow. 
The  maximum  value  of  the  relative  axial  velocity  does  not  follow  a  monotonic  trend  with 
L  because  as  the  length  becomes  small  the  precise  nature  of  the  beam  dynamics  vary 
in  a  complicated  manner  which  directly  affects  the  fluid  motion  and  therefore  the  fluid 
dissipation.  In  fact,  the  smallest  beam  L'/37.5  has  an  aspect  ratio  of  L/w  =  1  and  is  better 
described  as  a  plate  undergoing  complicated  dynamics  as  indicated  by  the  presence  of 
higher  mode  effects  in  Fig.  16(a).  Overall,  our  results  suggest  that  the  relative  magnitude 
of  the  axial  velocity  can  be  used  to  indicate  the  applicability  of  the  two-dimensional  theory. 
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Figure  17:  (color)  The  variation  of  the  quality  Q,  panel  (a),  and  of  the  resonant  frequency 
in  fluid  Ufluo,  panel  (b),  by  decreasing  the  length  L  of  the  beam  relative  to  the  value  of 
case  (1)  given  by  L'.  The  triangles  are  the  results  from  numerical  simulation  and  the  solid 
line  is  the  analytical  prediction. 

In  many  microscale  technologies  the  ability  to  sense  small  forces  is  important  and  there¬ 
fore  a  small  spring  constant  is  desirable.  In  light  of  this,  the  improved  performance,  as 
measured  by  increased  values  of  Q  and  ojf  with  increasing  w,  increasing  h  or  decreasing 
L  all  come  at  the  price  of  reduced  force  sensitivity.  Using  Eq.  (28)  to  estimate  k  yields 
its  dependence  upon  geometry  and  the  magnitude  of  the  improved  performance  follows 
the  same  trend  as  increasing  k.  Overall,  these  tradeoffs  would  need  to  be  balanced  in  a 
particular  application. 

The  stochastic  dynamics  of  micron  and  nanoscale  elastic  beams  can  be  directly  quantified 
using  deterministic  numerical  computations  for  the  precise  geometries  and  conditions  of 
experiment.  We  have  shown  that  the  geometry  of  doubly-clamped  beams  can  be  tailored 
to  overcome  the  strong  fluid  damping  that  occurs  for  small  scale  systems  in  a  viscous  fluid. 
Our  numerical  exploration  has  been  used  to  build  physical  insights  into  the  stochastic  dy¬ 
namics  and  to  place  realistic  bounds  upon  the  applicability  of  the  two-dimensional  theory. 
Overall,  we  find  that  the  two  dimensional  theory  is  quite  accurate  far  beyond  what  may 
have  been  expected  based  upon  the  underlying  assumptions.  When  deviations  do  occur  a 
significant  factor  are  fluid  velocities  in  the  axial  direction  resulting  in  increased  dissipa¬ 
tion  and  a  lower  added  mass.  It  is  anticipated  that  these  results  will  be  useful  in  guiding 
the  development  of  future  experiments  by  providing  the  basis  for  predictions  that  cover  a 
wide  range  of  geometries.  Furthermore,  our  results  provide  insight  into  the  development 
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Figure  18:  The  transverse  fluid  velocity  for  cases  1-8  along  the  line  beginning  at 
(0.0,0.0,0.01/rm)  and  ending  at  (L,0.0,0.01/im)  from  the  deterministic  numerical  simu¬ 
lations  where  the  beam  velocity  is  at  its  maximum  value.  The  baseline  geometry  is  shown 
as  the  dashed  line,  cases  2-7  are  the  solid  lines,  and  case  8  is  the  dash-dot  line. 

of  accurate  theoretical  models  valid  for  the  finite  geometries  used  in  experiment. 


5  Executive  Summary 

We  have  developed  analytical  and  theoretical  techniques  to  quantify  the  stochastic  and 
externally  driven  dynamics  of  elastic  objects  in  a  viscous  fluid  for  the  precise  conditions  of 
experiment.  These  techniques  allow  the  quantification  of  future  designs  aiming  to  exploit 
the  dynamics  of  micro  and  nanoscale  devices  in  fluid. 

This  project  provided  support  for  the  following  graduate  research.  The  theses  are  available 
to  the  public  in  digital  form  from  the  University  Library  at  Virginia  Tech  (www.lib.vt.edu). 
The  specific  thesis  details  are  given  below: 

•  Margarita  Villa,  Masters  Thesis,  Tailoring  the  geometry  of  micron  scale  resonators 
to  overcome  viscous  damping,  Virginia  Tech,  (2009). 

•  Matthew  Clark,  PhD  Dissertation,  The  driven  and  stochastic  dynamics  of  micro 
and  nanoscale  cantilevers  in  viscous  fluid  and  near  a  solid  boundary,  Virginia  Tech, 
(2008). 
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Figure  19:  The  normalized  axial  velocity  from  numerical  simulation,  (a)  The  axial  ve¬ 
locity  found  when  varying  the  width  and  height  of  the  baseline  geometry,  (b)  The  axial 
velocity  found  when  decreasing  the  length  of  the  baseline  geometry.  The  results  for  the 
baseline  geometry  are  given  by  the  dashed  line.  The  axial  velocities  are  computed  along 
a  line  in  the  x-direction  with  origin  (0,0,0.01/rm)  and  end-point  (L,0,0.01/im),  see  Fig.  1 1 
for  definitions  of  coordinate  directions.  The  velocity  is  normalized  by  mq  for  each  case 
where  uq  is  the  maximum  transverse  velocity  which  occurs  at  a;  =  L/2.  The  abscissa  is 
normalized  by  the  length  L  for  each  case. 

•  Carlos  Carvajal,  Masters  Thesis,  The  fluid  coupled  motion  of  micro  and  nanoscale 
cantilevers,  Virginia  Tech,  (2007). 
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